########################################
########################################
########################################
########################################
##Working directory should be set to location of this R file.
########################################
########################################
########################################
########################################

########################################
########################################
########################################
########################################
##Replication code for statistical analyses reported in "Symbols of State: Explaining Prestige Projects in the Global South"
##International Studies Quarterly
##Austin Strange
##Last updated: 6 March 2024
########################################
########################################
########################################
########################################

#Code for installing Zelig if archived and otherwise unavailable
#install.packages("https://cran.r-project.org/src/contrib/Archive/Zelig/Zelig_5.1.6.1.tar.gz", 
#      repos=NULL, 
#      type="source")

##Load packages
library(AER)
library(arsenal)
library(boot)
library(boot.pval)
library(brglm2)
library(broom)
library(data.table)
library(dplyr)
library(geepack)
library(ggplot2)
library(ggpubr)
library(gridExtra)
library(lmtest)
library(MCMCpack)
library(psych)
library(sandwich)
library(stargazer)
library(tidyr)
library(VGAM)
library(xtable)
library(Zelig)

########################################
########################################
########################################
########################################
##Load, merge, and clean data
########################################
########################################
########################################
######################################## 

d <- read.csv("ISQ_AMS_AllocationDataset.csv")

########################################
########################################
########################################
########################################
########################################
########################################
##Manuscript statistical figures and tables 
########################################
########################################
########################################
########################################
########################################
########################################

#Create subsetted data for additional study periods
d2 <- d[d$Year > 1999,]
d71 <- d[d$Year > 1970,]

########################################
########################################
#Table 1: Cross-national Allocation of Prestige Projects
########################################
########################################

##Table 1 (a) Small States and Prestige Project Allocation
t1m1 <- glm(P1B ~ L1PennGDP_2010usd_ln + factor(Year), data=d,family = "binomial", 
            method = "brglmFit"); summary(t1m1)
t1m2 <- glm(P1B ~ L1PennGDP_2010usd_ln + L1ROC + L1Democracy + english + factor(Year), data=d,family = "binomial", 
            method = "brglmFit"); summary(t1m2)
t1m3 <- glm(P1B ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
             L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d71,family = "binomial", 
            method = "brglmFit"); summary(t1m3)
t1m4 <- glm(P1B ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
              L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d2,family = "binomial", 
            method = "brglmFit"); summary(t1m4)
t1m5 <- lm(AP1_2010usd_ln ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
             L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year),data=d2); summary(t1m5)

cse_t1m1 <- coeftest(t1m1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t1m2 <- coeftest(t1m2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t1m3 <- coeftest(t1m3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t1m4 <- coeftest(t1m4,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t1m5 <- coeftest(t1m5,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_t1m1,cse_t1m2,cse_t1m3,cse_t1m4,cse_t1m5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = "Year",
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")))

##Table 1 (b) Political Support and Prestige Project Allocation
t1bm1 <- glm(P1B ~ L1daoge + factor(Host) + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(t1bm1)
t1bm2 <- glm(P1B ~ L1daoge + L1PennGDP_2010usd_ln + L1Democracy + english + factor(Host) + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(t1bm2)
t1bm3 <- glm(P1B ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
              L1unsc + L1dis.tot.aff_ln + english + L1toda_2010usd_ln + factor(Host) + factor(Year), data=d71,family = "binomial", 
             method = "brglmFit"); summary(t1bm3)
t1bm4 <- glm(P1B ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
              L1unsc + L1dis.tot.aff_ln + english + L1toda_2010usd_ln + factor(Host) + factor(Year), data=d2,family = "binomial", 
             method = "brglmFit"); summary(t1bm4)
t1bm5 <- lm(AP1_2010usd_ln ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
             L1unsc + L1dis.tot.aff_ln + english + L1toda_2010usd_ln + factor(Host) + factor(Year),data=d2); summary(t1bm5)

cse_t1bm1 <- coeftest(t1bm1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t1bm2 <- coeftest(t1bm2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t1bm3 <- coeftest(t1bm3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t1bm4 <- coeftest(t1bm4,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t1bm5 <- coeftest(t1bm5,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_t1bm1,cse_t1bm2,cse_t1bm3,cse_t1bm4,cse_t1bm5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = c("Year","Host"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")))

########################################
########################################
#Table 2: Comparing Prestige and Non-Prestige Project Allocation
########################################
########################################

##Table 2 Columns 1-3: Small States
t2am1 <- lm(AP1_2010usd_ln ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
               L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year),data=d2); summary(t2am1)
t2am2 <- lm(ANONPODA_2010usd_ln ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
              L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year),data=d2); summary(t2am2)
t2am3 <- lm(ANONPOOF_2010usd_ln ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
              L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year),data=d2); summary(t2am3)

cse_t2am1 <- coeftest(t2am1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t2am2 <- coeftest(t2am2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t2am3 <- coeftest(t2am3,vcov = vcovCL,type = "HC1",cluster = ~Host)

##Table 2 Columns 4-6: Political Support
t2bm1 <- lm(AP1_2010usd_ln ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy  +
              L1unsc + L1dis.tot.aff_ln + english + L1toda_2010usd_ln + factor(Host) + factor(Year),data=d2); summary(t2bm1)
t2bm2 <- lm(ANONPODA_2010usd_ln ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy +
             L1unsc + L1dis.tot.aff_ln + english + L1toda_2010usd_ln + factor(Host) + factor(Year),data=d2); summary(t2bm2)
t2bm3 <- lm(ANONPOOF_2010usd_ln ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy  +
             L1unsc + L1dis.tot.aff_ln + english + L1toda_2010usd_ln + factor(Host) + factor(Year),data=d2); summary(t2bm3)

cse_t2bm1 <- coeftest(t2bm1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t2bm2 <- coeftest(t2bm2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t2bm3 <- coeftest(t2bm3,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_t2am1,cse_t2am2,cse_t2am3,cse_t2bm1,cse_t2bm2,cse_t2bm3,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = c("Year","Host"),
          omit.stat=c("ser","f","adj.rsq"))

########################################
########################################
#Figure 3: Perceptions of Project Features for Prestige vs. Non-Prestige Projects
########################################
########################################

#Read in PNG survey data
png <- read.csv("ISQ_AMS_PNGSurveyDataset.csv")
png <- png[,c(2:14)]

pvis <- t.test(png$vis[png$p == 1], png$vis[png$p == 0], data = png) 
pvis_e <- as.data.frame(rbind(pvis$estimate))
pvis_ci <- as.data.frame(rbind(pvis$conf.int))
pvis2 <- cbind(pvis_e,pvis_ci); colnames(pvis2) <- c("mx","my","cil","cih"); pvis2$dm <- pvis2$mx-pvis2$my; pvis2$test <- "Visibility"; pvis2$type <- "PNG"

pngng <- t.test(png$attribute_natgov[png$p == 1], png$attribute_natgov[png$p == 0], data = png) 
pngng_e <- as.data.frame(rbind(pngng$estimate))
pngng_ci <- as.data.frame(rbind(pngng$conf.int))
pngng2 <- cbind(pngng_e,pngng_ci); colnames(pngng2) <- c("mx","my","cil","cih"); pngng2$dm <- pngng2$mx-pngng2$my; pngng2$test <- "Attribution"; pngng2$type <- "PNG"

pMot_IntStat <- t.test(png$Mot_IntStat[png$p == 1], png$Mot_IntStat[png$p == 0], data = png) 
pMot_IntStat_e <- as.data.frame(rbind(pMot_IntStat$estimate))
pMot_IntStat_ci <- as.data.frame(rbind(pMot_IntStat$conf.int))
pMot_IntStat2 <- cbind(pMot_IntStat_e,pMot_IntStat_ci); colnames(pMot_IntStat2) <- c("mx","my","cil","cih"); pMot_IntStat2$dm <- pMot_IntStat2$mx-pMot_IntStat2$my; pMot_IntStat2$test <- "Status"; pMot_IntStat2$type <- "PNG"

plot1a <- rbind(pvis2,pngng2,pMot_IntStat2)
plot1a$test <- factor(plot1a$test,levels=rev(unique(plot1a$test)))
  
p1 <- ggplot(plot1a, aes(y=dm, x=test)) +
  geom_point(position = position_dodge(width = .3),size=3,color="darkblue") +
  geom_errorbar(aes(ymax=cih, ymin=cil), color="dark gray", position = position_dodge(0.3), width=.1) +
  geom_hline(yintercept=c(0), linetype="dotted",size=1) +
  xlab("Project Feature") + ylab("Difference in Means") +
  theme_minimal() + 
  theme(axis.title.x = element_text(size=20),
        axis.text.x  = element_text(size=14, color="black"),
        axis.text.y  = element_text(size=14, color="black"),
        axis.title.y = element_text(size=20),
        legend.title=element_text(size=20),
        legend.text=element_text(size=14)) +
  labs(color = "Country") +
  guides(fill = "none",color = guide_legend(reverse=TRUE)) +
  coord_flip()

########################################
########################################
########################################
########################################
########################################
########################################
##Appendices with statistical figures and tables
########################################
########################################
########################################
########################################
########################################
########################################

########################################
########################################
########################################
########################################
##Appendix B: Additional Tests of Project Allocation
########################################
########################################
########################################
########################################

########################################
########################################
#Table A4: Country and Year Coverage in Statistical Analyses
########################################
########################################
earliest <- d %>% group_by(Host) %>% slice(which.min(Year)); earliest <- earliest[,c(4,2)]
latest <- d %>% group_by(Host) %>% slice(which.max(Year)); latest <- latest[,c(4,2)]; colnames(latest)[2] <- "Year2"
yrs <- merge(earliest,latest,by=c("Host"))
cols <- c("Year","Year2")
yrs$Range <- apply(yrs[ ,cols ],1,paste,collapse = "-")
yrs <- yrs[,c(1,4)]
yrs$Range <- paste0("(", yrs$Range, ")")
cols2 <- c("Host","Range")
yrs$Country <- apply(yrs[ ,cols2],1,paste,collapse = " ")
yrs <- as.data.frame(yrs[,c(3)])
yrs1 <- as.data.frame(yrs[c(1:42),])
yrs2 <- as.data.frame(yrs[c(43:84),])
yrs3 <- as.data.frame(yrs[c(85:125),])
yrs3[42,] <- "NA"
yrsall <- cbind(yrs1,yrs2,yrs3)
print(xtable(yrsall),include.rownames=FALSE)

########################################
########################################
#Table A5: Summary Statistics for Variables in Main Analyses
########################################
########################################
decstats <- describe(d[,c(5,8,11,12,15:27)],skew=F,na.rm=T)
decstats <- setDT(decstats, keep.rownames = TRUE)[]
xtable(decstats)

########################################
########################################
#Table A6: Cross-national Allocation of Prestige Projects (Maintenance Excluded)
########################################
########################################

##(a) Small States and Prestige Project Allocation
at1m1 <- glm(P12B ~ L1PennGDP_2010usd_ln + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(at1m1)
at1m2 <- glm(P12B ~ L1PennGDP_2010usd_ln + L1ROC + L1Democracy + english + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(at1m2)
at1m3 <- glm(P12B ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
              L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d71,family = "binomial", 
             method = "brglmFit"); summary(at1m3)
at1m4 <- glm(P12B ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
              L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d2,family = "binomial", 
             method = "brglmFit"); summary(at1m4)
at1m5 <- lm(AP12_2010usd_ln ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
             L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year),data=d2); summary(at1m5)

cse_at1m1 <- coeftest(at1m1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_at1m2 <- coeftest(at1m2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_at1m3 <- coeftest(at1m3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_at1m4 <- coeftest(at1m4,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_at1m5 <- coeftest(at1m5,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_at1m1,cse_at1m2,cse_at1m3,cse_at1m4,cse_at1m5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = "Year",
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")))

##(b) Political Support and Prestige Project Allocation
at2m1 <- glm(P12B ~ L1daoge + factor(Host) + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(at2m1)
at2m2 <- glm(P12B ~ L1daoge + L1PennGDP_2010usd_ln + L1Democracy + english + factor(Host) + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(at2m2)
at2m3 <- glm(P12B ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
              L1unsc + L1dis.tot.aff_ln + english + L1toda_2010usd_ln + factor(Host) + factor(Year), data=d71,family = "binomial", 
             method = "brglmFit"); summary(at2m3)
at2m4 <- glm(P12B ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
              L1unsc + L1dis.tot.aff_ln + english + L1toda_2010usd_ln + factor(Host) + factor(Year), data=d2,family = "binomial", 
             method = "brglmFit"); summary(at2m4)
at2m5 <- lm(AP12_2010usd_ln ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
             L1unsc + L1dis.tot.aff_ln + english + L1toda_2010usd_ln + factor(Host) + factor(Year),data=d2); summary(at2m5)
                 
cse_at2m1 <- coeftest(at2m1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_at2m2 <- coeftest(at2m2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_at2m3 <- coeftest(at2m3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_at2m4 <- coeftest(at2m4,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_at2m5 <- coeftest(at2m5,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_at2m1,cse_at2m2,cse_at2m3,cse_at2m4,cse_at2m5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = c("Year","Host"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")))

########################################
########################################
#Table A7: Cross-national Allocation of Prestige Projects (Maintenance and Renovation Excluded)
########################################
########################################

##(a) Small States and Prestige Project Allocation
rt1m1 <- glm(P13B ~ L1PennGDP_2010usd_ln + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(rt1m1)
rt1m2 <- glm(P13B ~ L1PennGDP_2010usd_ln + L1ROC + L1Democracy + english + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(rt1m2)
rt1m3 <- glm(P13B ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
               L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d71,family = "binomial", 
             method = "brglmFit"); summary(rt1m3)
rt1m4 <- glm(P13B ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
               L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d2,family = "binomial", 
             method = "brglmFit"); summary(rt1m4)
rt1m5 <- lm(AP13_2010usd_ln ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
              L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year),data=d2); summary(rt1m5)

cse_rt1m1 <- coeftest(rt1m1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_rt1m2 <- coeftest(rt1m2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_rt1m3 <- coeftest(rt1m3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_rt1m4 <- coeftest(rt1m4,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_rt1m5 <- coeftest(rt1m5,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_rt1m1,cse_rt1m2,cse_rt1m3,cse_rt1m4,cse_rt1m5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = c("Year","Host"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")))

##(b) Political Support and Prestige Project Allocation
rt2m1 <- glm(P13B ~ L1daoge + factor(Host) + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(rt2m1)
rt2m2 <- glm(P13B ~ L1daoge + L1PennGDP_2010usd_ln + L1Democracy + english + factor(Host) + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(rt2m2)
rt2m3 <- glm(P13B ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
               L1unsc + L1toda_2010usd_ln + L1dis.tot.aff_ln + english + factor(Host) + factor(Year), data=d71,family = "binomial", 
             method = "brglmFit"); summary(rt2m3)
rt2m4 <- glm(P13B ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
               L1unsc + L1toda_2010usd_ln + L1dis.tot.aff_ln + english + factor(Host) + factor(Year), data=d2,family = "binomial", 
             method = "brglmFit")
rt2m5 <- lm(AP13_2010usd_ln ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
              L1unsc +  L1toda_2010usd_ln + L1dis.tot.aff_ln + english + factor(Host) + factor(Year),data=d2); summary(rt2m5)

cse_rt2m1 <- coeftest(rt2m1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_rt2m2 <- coeftest(rt2m2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_rt2m3 <- coeftest(rt2m3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_rt2m4 <- coeftest(rt2m4,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_rt2m5 <- coeftest(rt2m5,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_rt2m1,cse_rt2m2,cse_rt2m3,cse_rt2m4,cse_rt2m5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = c("Year","Host"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")))

########################################
########################################
#Table A8: Cross-national Allocation of Prestige Projects (Composite Measure)
########################################
########################################

##(a) Small States and Prestige Project Allocation
t1GIm1 <- glm(GIB ~ L1PennGDP_2010usd_ln + factor(Year), data=d,family = "binomial", 
              method = "brglmFit"); summary(t1GIm1)
t1GIm2 <- glm(GIB ~ L1PennGDP_2010usd_ln + L1ROC + L1Democracy + english + factor(Year), data=d,family = "binomial", 
              method = "brglmFit"); summary(t1GIm2)
t1GIm3 <- glm(GIB ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
                L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d71,family = "binomial", 
              method = "brglmFit"); summary(t1GIm3)
t1GIm4 <- glm(GIB ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
                L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d2,family = "binomial", 
              method = "brglmFit"); summary(t1GIm4)
t1GIm5 <- lm(AGI_2010usd_ln ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
               L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year),data=d2); summary(t1GIm5)

cse_t1GIm1 <- coeftest(t1GIm1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t1GIm2 <- coeftest(t1GIm2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t1GIm3 <- coeftest(t1GIm3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t1GIm4 <- coeftest(t1GIm4,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t1GIm5 <- coeftest(t1GIm5,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_t1GIm1,cse_t1GIm2,cse_t1GIm3,cse_t1GIm4,cse_t1GIm5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = "Year",
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")))

####(b) Political Support and Prestige Project Allocation
t2GIm1 <- glm(GIB ~ L1daoge + factor(Host) + factor(Year), data=d,family = "binomial", 
              method = "brglmFit"); summary(t2GIm1)
t2GIm2 <- glm(GIB ~ L1daoge + L1PennGDP_2010usd_ln + L1Democracy + english + factor(Host) + factor(Year), data=d,family = "binomial", 
              method = "brglmFit"); summary(t2GIm2)
t2GIm3 <- glm(GIB ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
                L1unsc + L1dis.tot.aff_ln + english + L1toda_2010usd_ln + factor(Host) + factor(Year), data=d71,family = "binomial", 
              method = "brglmFit"); summary(t2GIm3)
t2GIm4 <- glm(GIB ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
                L1unsc + L1dis.tot.aff_ln + english + L1toda_2010usd_ln + factor(Host) + factor(Year), data=d2,family = "binomial", 
              method = "brglmFit"); summary(t2GIm4)
t2GIm5 <- lm(AGI_2010usd_ln ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
               L1unsc + L1dis.tot.aff_ln + english + L1toda_2010usd_ln + factor(Host) + factor(Year),data=d2); summary(t2GIm5)

cse_t2GIm1 <- coeftest(t2GIm1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t2GIm2 <- coeftest(t2GIm2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t2GIm3 <- coeftest(t2GIm3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t2GIm4 <- coeftest(t2GIm4,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_t2GIm5 <- coeftest(t2GIm5,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_t2GIm1,cse_t2GIm2,cse_t2GIm3,cse_t2GIm4,cse_t2GIm5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = c("Year","Host"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")))

########################################
########################################
#Table A9: Alternative Small State Measures: GDP Per Capita and Population
########################################
########################################
ss1m1 <- glm(P1B ~ L1PennGDPPC_2010usd_ln + L1PennPop_ln + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(ss1m1)
ss1m2 <- glm(P1B ~ L1PennGDPPC_2010usd_ln + L1PennPop_ln + L1ROC + L1Democracy + english + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(ss1m2)
ss1m3 <- glm(P1B ~ L1PennGDPPC_2010usd_ln + L1PennPop_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
              L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d71,family = "binomial", 
             method = "brglmFit"); summary(ss1m3)
ss1m4 <- glm(P1B ~ L1PennGDPPC_2010usd_ln + L1PennPop_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
              L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d2,family = "binomial", 
             method = "brglmFit"); summary(ss1m4)
ss1m5 <- lm(AP1_2010usd_ln ~ L1PennGDPPC_2010usd_ln + L1PennPop_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
             L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year),data=d2); summary(ss1m5)

cse_ss1m1 <- coeftest(ss1m1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_ss1m2 <- coeftest(ss1m2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_ss1m3 <- coeftest(ss1m3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_ss1m4 <- coeftest(ss1m4,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_ss1m5 <- coeftest(ss1m5,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_ss1m1,cse_ss1m2,cse_ss1m3,cse_ss1m4,cse_ss1m5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = c("Year","Host"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")))

########################################
########################################
#Table A10: Alternative Small State Measures: State Strength
########################################
########################################
ss2m1 <- glm(P1B ~ L1Strength2_2010usd_ln + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(ss2m1)
ss2m2 <- glm(P1B ~ L1Strength2_2010usd_ln + L1ROC + L1Democracy + english + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(ss2m2)
ss2m3 <- glm(P1B ~ L1Strength2_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
               L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d71,family = "binomial", 
             method = "brglmFit"); summary(ss2m3)
ss2m4 <- glm(P1B ~ L1Strength2_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
               L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d2,family = "binomial", 
             method = "brglmFit"); summary(ss2m4)
ss2m5 <- lm(AP1_2010usd_ln ~ L1Strength2_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
              L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year),data=d2); summary(ss2m5)

cse_ss2m1 <- coeftest(ss2m1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_ss2m2 <- coeftest(ss2m2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_ss2m3 <- coeftest(ss2m3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_ss2m4 <- coeftest(ss2m4,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_ss2m5 <- coeftest(ss2m5,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_ss2m1,cse_ss2m2,cse_ss2m3,cse_ss2m4,cse_ss2m5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = c("Year","Host"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")))

########################################
########################################
#Table A11: Alternative Small State Measures: Central Government Budgets
########################################
########################################
ss3m1 <- glm(P1B ~ L1SPEED_ln + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(ss3m1)
ss3m2 <- glm(P1B ~ L1SPEED_ln + L1ROC + L1Democracy + english + factor(Year), data=d,family = "binomial", 
             method = "brglmFit"); summary(ss3m2)
ss3m3 <- glm(P1B ~ L1SPEED_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
               L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d71,family = "binomial", 
             method = "brglmFit"); summary(ss3m3)
ss3m4 <- glm(P1B ~ L1SPEED_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
               L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d2,family = "binomial", 
             method = "brglmFit"); summary(ss3m4)
ss3m5 <- lm(AP1_2010usd_ln ~ L1SPEED_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
              L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year),data=d2); summary(ss3m5)

cse_ss3m1 <- coeftest(ss3m1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_ss3m2 <- coeftest(ss3m2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_ss3m3 <- coeftest(ss3m3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_ss3m4 <- coeftest(ss3m4,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_ss3m5 <- coeftest(ss3m5,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_ss3m1,cse_ss3m2,cse_ss3m3,cse_ss3m4,cse_ss3m5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = c("Year","Host"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")))

########################################
########################################
#Table A12: Cross-national Prestige Project Allocation (Rare Events Logit)
########################################
########################################
set.seed(93331)
rer <- zelig(P1B ~ L1PennGDP_2010usd_ln + L1ROC + L1Democracy + english + factor(Year),
              data = d, tau = 416/6685, model = "relogit", bootstrap=500); summary(rer)

########################################
########################################
#Table A13: Cross-national Allocation of Prestige Projects (2-Year Lags)
########################################
########################################
L2m1 <- glm(P1B ~ L2PennGDP_2010usd_ln + factor(Year), data=d,family = "binomial", 
            method = "brglmFit"); summary(L2m1)
L2m2 <- glm(P1B ~ L2PennGDP_2010usd_ln + L2ROC + L2Democracy + english + factor(Year), data=d,family = "binomial", 
            method = "brglmFit"); summary(L2m2)
L2m3 <- glm(P1B ~ L2PennGDP_2010usd_ln + L2total_trade_2010usd_ln + L2ResDep.pc + L2ROC + L2Democracy + 
              L2absidealdiff + L2unsc + english + L2dis.tot.aff_ln + L2toda_2010usd_ln + factor(Year), data=d71,family = "binomial", 
            method = "brglmFit"); summary(L2m3)
L2m4 <- glm(P1B ~ L2PennGDP_2010usd_ln + L2total_trade_2010usd_ln + L2ResDep.pc + L2ROC + L2Democracy + 
              L2absidealdiff + L2unsc + english + L2dis.tot.aff_ln + L2toda_2010usd_ln + factor(Year), data=d2,family = "binomial", 
            method = "brglmFit"); summary(L2m4)
L2m5 <- lm(AP1_2010usd_ln ~ L2PennGDP_2010usd_ln + L2total_trade_2010usd_ln + L2ResDep.pc + L2ROC + L2Democracy + 
             L2absidealdiff + L2unsc + english + L2dis.tot.aff_ln + L2toda_2010usd_ln + factor(Year),data=d2); summary(L2m5)

cse_L2m1 <- coeftest(L2m1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_L2m2 <- coeftest(L2m2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_L2m3 <- coeftest(L2m3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_L2m4 <- coeftest(L2m4,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_L2m5 <- coeftest(L2m5,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_L2m1,cse_L2m2,cse_L2m3,cse_L2m4,cse_L2m5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = c("Year","Host"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")))

########################################
########################################
#Table A14: Cross-national Allocation of Prestige Projects (3-Year Lags)
########################################
########################################
L3m1 <- glm(P1B ~ L3PennGDP_2010usd_ln + factor(Year), data=d,family = "binomial", 
            method = "brglmFit"); summary(L3m1)
L3m2 <- glm(P1B ~ L3PennGDP_2010usd_ln + L3ROC + L3Democracy + english + factor(Year), data=d,family = "binomial", 
            method = "brglmFit"); summary(L3m2)
L3m3 <- glm(P1B ~ L3PennGDP_2010usd_ln + L3total_trade_2010usd_ln + L3ResDep.pc + L3ROC + L3Democracy + 
              L3absidealdiff + L3unsc + english + L3dis.tot.aff_ln + L3toda_2010usd_ln + factor(Year), data=d71,family = "binomial", 
            method = "brglmFit"); summary(L3m3)
L3m4 <- glm(P1B ~ L3PennGDP_2010usd_ln + L3total_trade_2010usd_ln + L3ResDep.pc + L3ROC + L3Democracy + 
              L3absidealdiff + L3unsc + english + L3dis.tot.aff_ln + L3toda_2010usd_ln + factor(Year), data=d2,family = "binomial", 
            method = "brglmFit"); summary(L3m4)
L3m5 <- lm(AP1_2010usd_ln ~ L3PennGDP_2010usd_ln + L3total_trade_2010usd_ln + L3ResDep.pc + L3ROC + L3Democracy + 
             L3absidealdiff + L3unsc + english + L3dis.tot.aff_ln + L3toda_2010usd_ln + factor(Year),data=d2); summary(L3m5)

cse_L3m1 <- coeftest(L3m1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_L3m2 <- coeftest(L3m2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_L3m3 <- coeftest(L3m3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_L3m4 <- coeftest(L3m4,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_L3m5 <- coeftest(L3m5,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_L3m1,cse_L3m2,cse_L3m3,cse_L3m4,cse_L3m5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = c("Year","Host"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes")))

########################################
########################################
#Table A15: Cross-national Allocation of Prestige Projects (Institutions/Clientelism)
########################################
########################################
polim1 <- glm(P1B ~ L1v2xnp_client + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC  + 
              L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d71,family = "binomial", 
            method = "brglmFit"); summary(polim1)
polim2 <- glm(P1B ~ L1polity2 + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC  + 
                L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d71,family = "binomial", 
              method = "brglmFit"); summary(polim2)
polim3 <- glm(P1B ~ L1libdem + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC  + 
                L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d71,family = "binomial", 
              method = "brglmFit"); summary(polim3)
polim4 <- glm(P1B ~ L1Democracy + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC  + 
                L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + L1toda_2010usd_ln + factor(Year), data=d71,family = "binomial", 
              method = "brglmFit"); summary(polim4)

cse_polim1 <- coeftest(polim1,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_polim2 <- coeftest(polim2,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_polim3 <- coeftest(polim3,vcov = vcovCL,type = "HC1",cluster = ~Host)
cse_polim4 <- coeftest(polim4,vcov = vcovCL,type = "HC1",cluster = ~Host)

stargazer(cse_polim1,cse_polim2,cse_polim3,cse_polim4,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          single.row = T,
          omit = c("Year","Host"),
          omit.stat=c("ser","f","adj.rsq"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes")))

########################################
########################################
#Figure B2: Simulations for Difference in Coefficients Reported in Table 2
########################################
########################################
set.seed(93331)

boot1 <- function(data=d2, indices) {
  d <- data[indices, ]
  t3m1 <- lm(AP1_2010usd_ln ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
               L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + factor(Year),data=d)
  t3m2 <- lm(ANONPODA_2010usd_ln ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
               L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + factor(Year),data=d)
  return(coef(t3m1)["L1PennGDP_2010usd_ln"] - coef(t3m2 )["L1PennGDP_2010usd_ln"])
}
boot1r <- boot(data=d2, statistic=boot1, R=1000)
boot.ci(boot1r,type=c("norm","basic","perc"))
boot.pval(boot1r)

boot2 <- function(data=d2, indices) {
  d <- data[indices, ]
  t3m1 <- lm(AP1_2010usd_ln ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
               L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + factor(Year),data=d)
  t3m3 <- lm(ANONPOOF_2010usd_ln ~ L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1ROC + L1Democracy + 
               L1absidealdiff + L1unsc + english + L1dis.tot.aff_ln + factor(Year),data=d)
  return(coef(t3m1)["L1PennGDP_2010usd_ln"] - coef(t3m3)["L1PennGDP_2010usd_ln"])
}
boot2r <- boot(data=d2, statistic=boot2, R=1000)
boot.ci(boot2r,type=c("norm","basic","perc"))
boot.pval(boot2r)

boot3 <- function(data=d2, indices) {
  d <- data[indices, ]
  t4m1 <- lm(AP1_2010usd_ln ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
               L1unsc + L1dis.tot.aff_ln + english + factor(Host) + factor(Year),data=d)
  t4m2 <- lm(ANONPODA_2010usd_ln ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
               L1unsc + L1dis.tot.aff_ln + english + factor(Host) + factor(Year),data=d)
  return(coef(t4m1)["L1daoge"] - coef(t4m2 )["L1daoge"])
}
boot3r <- boot(data=d2, statistic=boot3, R=1000)
boot.ci(boot3r,type=c("norm","basic","perc"))
boot.pval(boot3r)

boot4 <- function(data=d2, indices) {
  d <- data[indices, ]
  t4m1 <- lm(AP1_2010usd_ln ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
               L1unsc + L1dis.tot.aff_ln + english + factor(Host) + factor(Year),data=d)
  t4m3 <- lm(ANONPOOF_2010usd_ln ~ L1daoge + L1absidealdiff + L1PennGDP_2010usd_ln + L1total_trade_2010usd_ln + L1ResDep.pc + L1Democracy + 
               L1unsc + L1dis.tot.aff_ln + english + factor(Host) + factor(Year),data=d)
  return(coef(t4m1)["L1daoge"] - coef(t4m3)["L1daoge"])
}
boot4r <- boot(data=d2, statistic=boot4, R=1000)
boot.ci(boot4r,type=c("norm","basic","perc"))
boot.pval(boot4r)

BootDist.graph1 <- data.frame(xbar=boot1r$t)
bs1 <- ggplot(BootDist.graph1, aes(x=xbar)) +
  geom_histogram(bins=40) +
  ggtitle('Difference in Coefficients, Table 2 Columns 1 and 2') + 
  labs(y="Count",x="Difference in coefficient") +
  theme(plot.title = element_text(size=18,face="bold"),
        axis.title.x = element_text(size=18),
        axis.title.y = element_text(size=18),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14)) +
  ylim(0,82)

BootDist.graph2 <- data.frame(xbar=boot2r$t)
bs2 <- ggplot(BootDist.graph2, aes(x=xbar)) +
  geom_histogram(bins=40) +
  ggtitle('Difference in Coefficients, Table 2 Columns 1 and 3') + 
  labs(y="Count",x="Difference in coefficient")+
  theme(plot.title = element_text(size=18,face="bold"),
        axis.title.x = element_text(size=18),
        axis.title.y = element_text(size=18),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14))+
  ylim(0,82)

BootDist.graph3 <- data.frame(xbar=boot3r$t)
bs3 <-ggplot(BootDist.graph3, aes(x=xbar)) +
  geom_histogram(bins=40) +
  ggtitle('Difference in Coefficients, Table 2 Columns 4 and 5') + 
  labs(y="Count",x="Difference in coefficient")+
  theme(plot.title = element_text(size=18,face="bold"),
        axis.title.x = element_text(size=18),
        axis.title.y = element_text(size=18),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14))+
  ylim(0,82)

BootDist.graph4 <- data.frame(xbar=boot4r$t)
bs4 <-ggplot(BootDist.graph4, aes(x=xbar)) +
  geom_histogram(bins=40) +
  ggtitle('Difference in Coefficients, Table 2 Columns 4 and 6')+ 
  labs(y="Count",x="Difference in coefficient")+
  theme(plot.title = element_text(size=18,face="bold"),
        axis.title.x = element_text(size=18),
        axis.title.y = element_text(size=18),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14))+
  ylim(0,82)

grid.arrange(bs1,bs2,bs3,bs4,ncol=2)

########################################
########################################
########################################
########################################
##Appendix D: Supplemental Information for Survey Experiment
########################################
########################################
########################################
########################################

########################################
########################################
#Table A16: OLS Regressions for PNG Prestige Project Attributes
########################################
########################################
r1 <- lm(vis ~ Treatment_DO + age+gender + income2 + highered + national + aidknow + pretreat_chn2,data=png); summary(r1)
r2 <- lm(png$attribute_natgov ~ Treatment_DO + age+gender + income2 + highered + national + aidknow + pretreat_chn2,data=png); summary(r2)
r3 <- lm(Mot_IntStat ~ Treatment_DO + age+gender + income2 + highered + national + aidknow + pretreat_chn2,data=png); summary(r3)

stargazer(r1,r2,r3,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          omit.stat=c("ser","f","adj.rsq"))

########################################
########################################
#Figure D3: Differences in Means by Public- vs. Private-facing Prestige Projects
########################################
########################################
png$t <- png$Treatment_DO
png$t[png$t == "nonp_1"] <- "np"
png$t[png$t == "nonp_2"] <- "np"
png$t <- as.factor(png$t); png$t <- relevel(png$t,ref="np")

pvisA <- t.test(png$vis[png$t == "p_pub"], png$vis[png$t == "np"], data = png) 
pvisA_e <- as.data.frame(rbind(pvisA$estimate))
pvisA_ci <- as.data.frame(rbind(pvisA$conf.int))
pvisA2 <- cbind(pvisA_e,pvisA_ci); colnames(pvisA2) <- c("mx","my","cil","cih"); pvisA2$dm <- pvisA2$mx-pvisA2$my; pvisA2$test <- "Visibility"; pvisA2$type <- "Theater"

pvisB <- t.test(png$vis[png$t == "p_priv"], png$vis[png$t == "np"], data = png) 
pvisB_e <- as.data.frame(rbind(pvisB$estimate))
pvisB_ci <- as.data.frame(rbind(pvisB$conf.int))
pvisB2 <- cbind(pvisB_e,pvisB_ci); colnames(pvisB2) <- c("mx","my","cil","cih"); pvisB2$dm <- pvisB2$mx-pvisB2$my; pvisB2$test <- "Visibility"; pvisB2$type <- "Ministry"

pngngA <- t.test(png$attribute_natgov[png$t == "p_pub"], png$attribute_natgov[png$t == "np"], data = png) 
pngngA_e <- as.data.frame(rbind(pngngA$estimate))
pngngA_ci <- as.data.frame(rbind(pngngA$conf.int))
pngngA2 <- cbind(pngngA_e,pngngA_ci); colnames(pngngA2) <- c("mx","my","cil","cih"); pngngA2$dm <- pngngA2$mx-pngngA2$my; pngngA2$test <- "Attribution"; pngngA2$type <- "Theater"

pngngB <- t.test(png$attribute_natgov[png$t == "p_priv"], png$attribute_natgov[png$t == "np"], data = png) 
pngngB_e <- as.data.frame(rbind(pngngB$estimate))
pngngB_ci <- as.data.frame(rbind(pngngB$conf.int))
pngngB2 <- cbind(pngngB_e,pngngB_ci); colnames(pngngB2) <- c("mx","my","cil","cih"); pngngB2$dm <- pngngB2$mx-pngngB2$my; pngngB2$test <- "Attribution"; pngngB2$type <- "Ministry"

pMot_IntStatA <- t.test(png$Mot_IntStat[png$t == "p_pub"], png$Mot_IntStat[png$t == "np"], data = png) 
pMot_IntStatA_e <- as.data.frame(rbind(pMot_IntStatA$estimate))
pMot_IntStatA_ci <- as.data.frame(rbind(pMot_IntStatA$conf.int))
pMot_IntStatA2 <- cbind(pMot_IntStatA_e,pMot_IntStatA_ci); colnames(pMot_IntStatA2) <- c("mx","my","cil","cih"); pMot_IntStatA2$dm <- pMot_IntStatA2$mx-pMot_IntStatA2$my; pMot_IntStatA2$test <- "Status"; pMot_IntStatA2$type <- "Theater"

pMot_IntStatB <- t.test(png$Mot_IntStat[png$t == "p_priv"], png$Mot_IntStat[png$t == "np"], data = png) 
pMot_IntStatB_e <- as.data.frame(rbind(pMot_IntStatB$estimate))
pMot_IntStatB_ci <- as.data.frame(rbind(pMot_IntStatB$conf.int))
pMot_IntStatB2 <- cbind(pMot_IntStatB_e,pMot_IntStatB_ci); colnames(pMot_IntStatB2) <- c("mx","my","cil","cih"); pMot_IntStatB2$dm <- pMot_IntStatB2$mx-pMot_IntStatB2$my; pMot_IntStatB2$test <- "Status"; pMot_IntStatB2$type <- "Ministry"

plot2 <- rbind(pvisA2,pvisB2,pngngA2,pngngB2,pMot_IntStatA2,pMot_IntStatB2)

plot2$test <- factor(plot2$test,levels=rev(unique(plot2$test)))
plot2$type <- as.factor(plot2$type)

p2 <- ggplot(plot2, aes(y=dm, x=test,fill=type,color=type)) +
  scale_fill_viridis_d() +
  geom_point(position = position_dodge(width = .3),size=3) +
  geom_errorbar(aes(ymax=cih, ymin=cil), color="dark gray", position = position_dodge(0.3), width=.1) +
  geom_hline(yintercept=c(0), linetype="dotted",size=1) +
  xlab("Project Feature") + ylab("Difference in Means") +
  theme_minimal() + 
  theme(axis.title.x = element_text(size=20),
        axis.text.x  = element_text(size=14, color="black"),
        axis.text.y  = element_text(size=14, color="black"),
        axis.title.y = element_text(size=20),
        legend.title=element_text(size=20),
        legend.text=element_text(size=14)) +
  labs(color = "Project") +
  guides(fill = "none",color = guide_legend(reverse=TRUE)) +
  coord_flip() 
